Long transient dynamics in the Anderson-Holstein model out of equilibrium 
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We calculate the time dependent nonequilibrium current through a single level quantum dot 
strongly coupled to a vibrational mode. The nonequilibrium real time dynamics caused by an 
instantaneous coupling of the leads to the quantum dot is discussed using an approximate method. 
The approach, which is specially designed for the strong polaronic regime, is based on the so-called 
polaron tunneling approximation. Considering different initial dot occupations, we show that a 
common steady state is reached after times much larger than the typical electron tunneling times 
due to a polaron blocking effect in the dot charge. A direct comparison is made with numerically 
exact data, showing good agreement for the time scales accessible by the diagrammatic Monte Carlo 
simulation method. 

PACS numbers: 73.63.-b, 71.38.-k, 72.15.Qm 



I. INTRODUCTION 



Experimental progress in the last few years has en- 
abled a detailed study of transport phenomena in single- 
molecule junctions. Such a junction can be considered 
as a quantum dot contacted to two electrodes via a tun- 
neling coupling. Applying a finite voltage the electrons 
tunnel through the quantum dot. The charging of the 
molecule leads to elastic deformations of its geometry 
which causes a coupling between electronic and vibra- 
tional degrees of freedom. This gives rise to effects like 
steps in the current-voltage characteristics^"— and the 
formation of sidebands in the excitation spectraiiii^rii 
In general such a quantum dot setup can be described by 
the Anderson-Holsteini^ii^ model. When one is mainly 
interested in the effects caused by the vibrational mode 
of the molecule it is customary to consider a linearly cou- 
pled local phonon mode and a single level quantum dot 
with spinless electrons. 

Depending on the temperature T, the coupling strength 
of the electrodes to the dot F, the electron-phonon in- 
teraction A as well as the phonon frequency ujq, different 
physical regimes can be distinguished. In the classical 
regime, T ^ F, the problem can be treated with semi- 
classical approaches using master equations^ii^. On the 
other hand, in the quantum regime, T <g; F, a great the- 
oretical effort has been made to develop methods to de- 
scribe transport phenomena within this model. This in- 
cludes different approximate approaches (see for example 
Refs. [8l fl9l - l28j and references therein) as well as numer- 
ically exact methods such as the diagrammatic Monte 
Carlo simulation (diagMC) 14,29- ^^ auxiliary-field quan- 
tum Monte CarloS^ or the multilayer multiconfiguration 
time dependent Hartree method^i^ as well as the itera- 
tive path integral summation scheme^. 
While most of these methods address the steady state 



behavior of the system, the way how this steady state is 
built up is not yet well understood. In the strong po- 
laronic regime the transient is only addressed by mean 
field studies'^^ or numerical methods. Recent numerical 
calculations'^^ of the transient dynamics of the current in 
the Anderson-Holstein model were indicating that in this 
model there exists a large time scale over which differ- 
ent initial preparations lead to different transport prop- 
erties which might even lead to a bistable situation. A 
bistable behavior was predicted previously for the steady 
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as well as the time dependentis^ regime within 
a mean field theory. 

In this paper we address the transient behavior of the 
Anderson-Holstein model for the strong polaronic regime 
at T = 0. Our aim is to develop an approximate method 
in order to analyze the behavior at very long times 
which is inaccessible by numerically exact approaches. 
We are specially interested in understanding how the 
steady state is built up from different initial conditions 
for the dot occupation. For this purpose we extend a re- 
cently proposed strong coupling approximation, namely 
the polaron tunneling approximation (PTA)^"^, to non- 
stationary situations. The results obtained from this ap- 
proximation are in very good agreement with the numer- 
ical ones obtained by diagMC for the times accessible to 
the exact method. On the other hand, our approxima- 
tion shows how the system converges into a steady state 
solution for much larger time scales, regardless of the ini- 
tial condition. 

The structure of the paper is the following: In section 
iniwe briefly introduce the Anderson-Holstein model and 
then discuss our approach to the time dependent prob- 
lem in the strongly polaronic regime in section IIIII The 
results are discussed in section HVl where we first address 
time scales accessible by numerical methods so that a 
comparison can be made. Finally the long time regime is 
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discussed, and we give a simple interpretation of the po- 
laron blocking mechanism. The paper is closed by some 
concluding remarks. 

II. THE MODEL 

The setup consists of two electrodes, left (L) and right 
(R), which are contacted by a tunneling coupling to an 
atom or molecule (hereafter called quantum dot) which is 
modeled by a single electronic level. This level is coupled 
to a single phonon mode which can be considered as the 
most relevant vibrational mode of the atom or molecule. 
Such a system can be described by the spinless Anderson- 
Holstein model given by^lii^ 

H = Hu+ Hi^R + i/ph + Ht + Hi, (1) 

where = ^ud^d describes the quantum dot energy 
level Ed, where and d are the electron creation and an- 
nihilation operators on the dot. iJLR = J2a,k ^ak°^ak'^ak 
corresponds to the non-interacting leads, with a = L,R 
denoting the left and right electrode respectively. The 
electronic creation and annihilation operators on elec- 
trode a at energy level eak are denoted by a}^^. and aak- 
The two electrodes are kept at different chemical poten- 
tials via a constant bias voltage eV — — Mr, inducing 
a nonequilibrium current through the dot. The tunneling 
coupling to the dot is described by 

i?T = ^7a(aLd + H.c.) , (2) 

(X,k 

where 7a are the tunneling amplitudes. We additionally 
introduce the tunneling rates Fq, = 'K^'^pa, where pa is 
the density of states of electrode a which is assumed to 
be independent of the energy. We further assume Fl = 
Fr = F/2. 

The phonon mode is described by i/ph = ojob^h^ which 
models a vibrational degree of freedom of the molecule 
with a normal mode of frequency wq- The coupling of 
this mode to the dot level is described by 

Hi = \dU{b + b^) , (3) 

with the coupling constant A. Throughout this paper we 
set h ^ e = rrif, = 1. 

III. GREEN FUNCTION APPROACH FOR THE 
POLARONIC REGIME 

Despite its simple structure, no exact analytical solu- 
tion of the Anderson-Holstein model is known for arbi- 
trary parameters. Only in the limits of either a vanishing 
dot- lead coupling^--, often called the atomic limit, or in 
the absence of phonons^ an exact solution can be ob- 
tained. For the stationary case an exact solution can 



also be found in the limits en/F ±cx3^. 
A common procedure to address the steady state is to 
assume a decoupled situation in the infinitely past, that 
is, at t — — oo. Therefore, for any time of interest, the 
system is in its steady state and transient effects do not 
need to be treated explicitly, which often simplifies the 
calculations. Despite the success of this procedure the 
information about how the steady state is established 
cannot be gathered. Such transients can be studied as- 
suming an initially vanishing tunneling coupling between 
the dot and the electrodes. Then the tunneling coupling 
of the dot to the leads is switched on at a certain ini- 
tial time t — Q and the subsequent time evolution to the 
steady state can be analyzed. In the present work we 
develop an analytical approach based on nonequilibrium 
Green's function techniques to study this transient be- 
havior of the Anderson-Holstein model. 
In order to access the nonequilibrium properties of the 
dot, the Keldysh Green's functions have to be determined 

D{t, t') = -i {rcd{t)dHt')X{t)X\t')) , (4) 

where 7c denotes the time ordering operator on the 
Keldysh contour C. The averaging is performed with 
respect to the complete quantum mechanical state of 
the system. X = ^Vsib'' -b) -j-j^g phonon cloud op- 
erator, which is obtained from a unitary Firsov-Lang 
transformation^^ where g = (A/wq)^. This parameter is 
a measure of the number of phonons forming the phonon 
cloud. 

Due to the internal symmetries in Keldysh space"*^ it 
is sufficient to consider only the advanced D^(t,t'), re- 
tarded D'^{t,t'), and lesser D^{t,t') dot's Green's func- 
tion to describe the transient current and dot occupation. 
A closed form solution of the complete Green's function 
in Eq. ^ is hard to achieve since all diagrams containing 
multi-phonon correlations have to be evaluated explicitly. 
Therefore, a simple approximation is desirable in order 
to describe strong electron-phonon couplings in a generic 
non-stationary regime. 

In this manuscript the time dependence is addressed by 
a diagrammatic expansion in terms of the dot-lead cou- 
pling. The average of the Green's function in Eq. (U) is 
then defined with respect to some given initial prepara- 
tion. For the initial preparation two possible dot occu- 
pations are considered: Either the dot is empty, so that 
nY){t = 0) = 0, or occupied, nY){t = Q) = 1. 



A. Atomic limit 

A good reference for analyzing the strong coupling 
regime is provided by the atomic limit, defined as the 
limit in which the tunneling rate F tends to zero. Fol- 
lowing lisj the model can be solved exactly so that e. g. 
the atomic retarded dot's Green's function at zero tem- 
perature is given by 
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(5) 



where the polaron shifted energy level of the dot is 
Ed = Ed — A^/wo and the dot occupation is denoted by 
noit) = {S{t)d{t)). Notice that for a strictly isolated dot 
the charge is constant and can only take the values or 1. 
But when considered as the limiting case of F 0, n-oit) 
corresponds to the mean dot occupation of the coupled 
system32.. This consideration is useful in the following 
discussion about the self-consistent determination of the 
dot charge. 

In frequency domain the retarded atomic Green's func- 
tion has the form 
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where ■& is an infinitesimal. 
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(6) 



B. Approximated PTA 
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FIG. 1: First three Feynman diagrams of the Dyson series 
for the PTA. The decoupled dot's Green's function Ddot{t,t') 
is represented with a solid line, the decoupled leads' Green's 
function gc,{t,t') with dashed lines. The two point phonon 
correlator connecting two time integration variables t,t' is 
denoted by Aa = {TcX{t)X'^ {t')) . 

A first step for going beyond the atomic limit is pro- 
vided by the PTA^ which is based on the assumption 
that the phonons are instantaneously excited once the 
electron tunnels to the dot and deexcited right after the 
electron leaves it. Basically, the time scale of an elec- 
tron on the dot is given by t^i oc r~^. The inverse of 
the energy due to the polaron formation determines the 
time scale for an (de)excitation of the polaron, that is, 
Tph oc (A^/wo)^^- In the polaronic regime, for A F 
and \/ijOq ^ 1 , the lifetime of the electron on the dot is 
much larger than the (de)excitation time of the polaron 
Tel 3> Tph, so that the PTA becomes a reasonable approx- 
imation. 



The corresponding Feynman diagrams of the PTA Dyson 
equation are plotted in Fig. [TJ Essentially, the PTA re- 
places multi-polaron correlations by a series of two point 
correlators A2 = (7cX{t)X''{t')'), that is, all phonon pro- 
cesses with more than one polaron are neglected. 
In the steady state situation the Dyson equation in the 
Keldysh matrix representation can be solved in frequency 
space 



DpTA — Dat + DatSDpTA , 



(7) 



by inserting the atomic limit Green's function and the 
leads' self-energy 



(8) 



where cr^ is a Pauli matrix in Keldysh space and ga de- 
notes the Green's functions of the decoupled leads. In 
this way, e.g. the retarded Green's function can be cal- 
culated as 
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l + iTDl,{uj) 



(9) 



From this expression it is straightforward to calculate the 
self-consistent dot charge using 

no = ^ ^ / dujUuj)lmD'{uj) , (10) 

where fa{'^) denotes the Fermi distribution on the lead a. 
The corresponding mean current is then obtained from 



(11) 



It is important to notice that the value of the steady- 
state current in general depends on the value of the mean 
charge ud . In a non-selfconsistent approach in which the 
dot charge is assumed to be either or 1 for calculat- 
ing the retarded Green's function in Eq. © one would 
obtain two different values for the stationary current, as 
illustrated in Fig. [2j However, the self-consistent calcu- 
lation yields a current- voltage characteristic which lies in 
between these two results, thus implying the absence of 
bistability within this approach. A very special situation 
is that of an electron- hole symmetric case ed = 0, with a 
symmetric voltage drop. Integrating over this symmetric 
voltage window in Eq. (|lip gives an / — characteris- 
tic which does not depend actually on the mean charge 
(shown as an inset in Fig. [2]). In the non-symmetric case 
not only the actual value of the current is different, but 
also such main features such as the height of the steps at 
multiples of the phonon frequency and the length of the 
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the cross terms arising from mixing different multiphonon 
resonances. 
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FIG. 2: Steady state current vs. applied bias voltage cal- 
culated within the PTA with nu ~ (brown dashed lines), 
riD = 1 (black dotted lines) and with the self-consistent charge 
(blue lines) determined by Eq. (|10|) . The dot is off-resonant 
with eu = —for, A = 16T, lvq = ST. Inset: The same plot 
but for £0=0. 



FIG. 3: Spectral density Ao{uj) of the PTA with fully treating 
the dot occupation (blue lines) vs. APTA (red lines) for en ~ 
-lOr, A = ler, ujq = Sr and V = lOr. inset: Zoom of 
the figure showing the difference between APTA and PTA 
spectral densities. 



plateaus between two phonon steps. 
In spite of the simplicity of the expression for the PTA 
Green's fmrction (Eq. its generalization to a time- 
dependent situation is rather involved. A further sim- 
plification of the approximation in the limit g ^ 1 and 
loq ^ T can be achieved by noticing that the polaronic 
(multi-phonon) side-bands essentially do not overlap (see 
Fig. [3]). This allows to approximate the poles in the re- 
tarded PTA Green's function as independent Lorentzian 
functions 



D 



APTA 
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(12) 



where, in order to fit the correct broadening around each 
resonance, the parameters P,^ have to have the form 



f ; = f ; [(1 - ^d) + nuSi^o] 
f,+ = f / [riD + (1 - nu)Si^o] 



with 
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(13) 
(14) 



(15) 



This approximated PTA (APTA) gives a much simpler 
form for the retarded Green's function while all its main 
features are still preserved. See Fig. [3] for a comparison of 
the spectral densities from APTA and PTA in the strong 
coupling regime. 

We further notice that in general the APTA fulfills the 
PTA Dyson equation (Eq. ([7])) approximately, ignoring 



C. Time dependent APTA 

In this subsection, we provide an ansatz for describing 
the time dependent nonequilibrium transient current for 
two different initial occupations within the spirit of the 
PTA approach. In principle, one has to solve the time- 
dependent Dyson equation 

DAPTA(i,i') = Dat(i,t') + J dsi J ds2 

Dat(i, si) S(si, S2) Dapta(s2, i') , (16) 

where the coupling of the dot to the leads is switched on 
abruptly at t = via j{t) — 9{t)^. The corresponding 
self energies then are given by 

S(i, t') = e{t)e{t')^^a, [gL(t, t') + gR(t, t')] a, . (17) 

A full self-consistent solution of these integral equations 
is a formidable task. The main idea of our ansatz is to 
perform a quasi-adiabatic approach in which the charge 
of the dot is assumed to evolve slowly in time while 
the spectral density adapts to the instantaneous value 
of this charge. This is consistent with the general PTA 
picture, since the average occupation of the quantum 
dot changes on a time scale given by P"^ whereas the 
phononic degrees of freedom adapt to the dot occupa- 
tion on a time scale given by (A^/wo)^^- This permits 
a reasonable "closed" or compact solution of the dy- 
namical problem. Technically, we shall neglect the tran- 
sient effects in the retarded Green's functions while fo- 
cusing on the transient properties of the lesser Green's 
function, which allows to determine the time-dependent 
charge self-consistently. More explicitly, for the retarded 
Green's function our ansatz for t, > is 
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1=0 



(18) 



r 



where, similarly to the steady state, the side-band broad- 
enings have the form 

f;-(t)=f/[(l-nD(t))+nD(t)(5;=o] , (19) 
f+(t) = ti [nD(t) + (1 - nu{t))Si=o] , (20) 

with nD{t) being the instantaneous mean charge which 
has to be determined self-consistently from the lesser 
Green's function 



"D(t) = -iL'APTA(*;^) 



(21) 



This Green's function satisfies the corresponding Dyson 
equation 



-^APTA ^ (1 + -DaftA^"^) (1 + ^^-^APTa) 



-Capta^^-^apta : 



(22) 



where integration over internal time arguments is im- 
plicitly assumed. The advanced dot's Green's function 



J 



needed for this Dyson equation can be obtained from 
the retarded one by the general relation D^{t,t') = 
{U^ {t' ,t))* . In Eq. the initial condition is provided 
by Dq which is determined by the initial dot occupation 



- t') = ic-»e-'''°(*-*')nD(0)e9°" 



jo(t-t') 



(23) 



Finally, the self-energies in Eq. (|22p are given by 



Y.''{t,t') ^ -i6{t)T6{t-t') , (24) 
S<(t,O=i^(0^(t')^E /dc.e-"(*-*')/a(w). (25) 



Next, the two possible initial occupations nD(0) — and 
"-D (0) — 1 are discussed separately. In the first case only 
the second term on the rhs of Eq. (|22p contributes and 
the resulting time dependent dot occupation is 



(t) = ^ ^ j duof^{u)\{l-nj,{t))S-{uo,t)+n^{t)S+{uo,t)\\ 

T TD ^ 



(26) 



a=L,R 



r 



where 
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-\{u)-fD)t 



-Tf{t)tp±\uoH 



w - ED ± (jjqI + ir,'*' (t) 



(27) 



On the other hand, for the case when nD(0) = 1, there is 
an extra contribution (5?id(0 arising from the first term 
on the rhs of Eq. ([^^ given by 

oo / 
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X \l + il~nuit))A+{t) + nuit)A-{t)\ 



(28) 



where 



°° -1 _ -f±(t)t -iwo(/Tni)t 

A± (t ) = i y f,„ ^ . (29) 

' ^ ^ ^ " (ZTm)a;o-ir±(0 ^ ' 



m=0 



Reaching a unique stationary state would require that 
6TiD{t) — i> for t -> oo. Although this would be war- 
ranted in an exact time-dependent self-consistent PTA, 
the approximations done within APTA yield a small fi- 
nite correction which vanishes as F/wq tends to zero. In 
order to numerically evaluate the dot occupation, a finite 
time step is chosen and the equations for the dot occupa- 
tions (Eqs. ((26)) and ([28)) ) are solved iteratively starting 
from the initial condition nj:,{t = 0) = {0, 1}. 
The average current is calculated using the relation 



oo 

= 7'Re J dsD'^it, s) [glis, t) - gg(s, t)] , (30) 



where g^{t,t') denotes the Keldysh Green's function of 
the decoupled lead a. Inserting the Green's functions 
and performing the time integration one obtains 
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e-r, (*)t _ cj) sin ((wq/ - cj) - fl{t) cos ((wo^ - w) t)^ 



-V/2-J:d 



+nD(0- 



-r,+ (t)t (^(^^i + sin ((wq/ + w) i) - f cos ((wq/ + w) t) 



> . 



(31) 



Note, that the case of a non-interacting electronic quan- 
tum dot can be obtained by setting g Q, with the 
same result for the current and the dot occupation as in 
reference nM. 



IV. RESULTS 

In this section the time dependent APTA results are 
discussed and compared with the numerically exact data 
from the diagrammatic Monte Carlo (diagMC) simula- 
tion method. This method uses a diagrammatic expan- 
sion in the dot-lead tunneling coupling. The occurring 
time integrals are evaluated stochastically using a Monte 
Carlo algorithm. Time dependent observables such as 
current or dot occupation can be calculated for arbi- 
trary system parameters like coupling strength, voltage 
and temperature. Details of the diagMC method can be 
found for example in [l3l30l |. 

Despite being numerically exact the diagMC has a draw- 
back since it suffers from the so-called "sign problem": 
The stochastic Monte Carlo sum has to be performed 
over terms with alternating signs. This leads to large sta- 
tistical errors in the observables causing the CPU time to 
scale exponentially with the system time. Therefore, for 
any realistic setup it is only possible to simulate the time 
dependent observables up to system times of the order of 

lor-i. 

In order to check the reliability of our approach in sub- 
section IIV Al a comparison is made between the APTA 
and diagMC for times which are accessible by the latter 
method. The long time scales are discussed in subsection 
llYBl 

In all cases we show results for both the initially empty 
and occupied dot. 



A. Short time scales 

The APTA is expected to be valid in the strong po- 
laronic regime with not too large applied voltages where 
many-polaron correlations should be small. Therefore, 
we perform the comparison between APTA and diagMC 
in the polaronic regime with A = IGF, ujq — SF and 
eD = — lOF. The choice of these parameters is also 
guided by the observations of Ref . [s^l suggesting a strong 
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FIG. 4: Comparison between the current from diagMC (sym- 
bols) and APTA (straight lines) for ?d = -lOF, A = ioF, 
ojQ = 8T, V = 2r. The current from the initially empty (oc- 
cupied) dot are highlighted in red (green) color for the APTA 
and represented by dots (diamonds) for the diagMC. Inset: 
Same plot but for larger times. 



"bistable" like behavior of the system for this case. The 
voltages are increased from small values to large ones 
where inelastic processes, not included in the PTA pic- 
ture, become important. 

The transient currents with V ~ 2T for the two differ- 
ent initial occupations of the dot, empty or occupied, are 
plotted in Fig. 21 where a remarkable agreement between 
the APTA and the diagMC is observed. The APTA de- 
scribes the main transient behavior: The current from 
APTA is forming plateaus with a constant current which 
are followed by short time intervals with a rapid change. 
These two situations exchange each other with a period 
only depending on the frequency of the phonon. These 
large oscillations of the current can be interpreted as a 
shake up process due to the sudden connection of the 
leads to the dot. For larger times the phonon cloud re- 
laxes and the oscillations become gradually smaller. 
The transient dynamics in Fig. 2] are quite different for 
the two possible initial configurations, empty or occu- 
pied. Depending on the initial configuration one observes 
a peak or a dip at t = where n is an integer. For 

times t > 6T~^ both transient currents oscillate around 
their joint steady state value (see inset of Fig.[S]). There- 
fore, the time scale for reaching a unique steady state is 



of the order of several T ^. 
0.4 
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FIG. 5: The currents from diagMC and APTA are compared 
with the same color code and parameters as in Fig.[3]but with 
different voltages V = 5T (figure on the top) and V = 26T 
(bottom). The steady state current cannot be reached even 
for times of the order of lOF"^. 

Increasing the voltage to V — dT the time dependent 
behavior changes: The transient current of the initially 
empty dot is much larger than the current for the ini- 
tially occupied one even for the largest times accessed 
by diagMC (see upper plot in Fig. [5|). Clearly, no joint 
steady state will be reached within times of the order of 
lOr^^ which is different to the observations for other 
quantum dot systems such as the Anderson impurity 
model*^. Therefore, this effect can be identified as a pure 
phononic one. The phonons in this regime seem to block 
the current depending on the initial configuration as it 
was shown previously in Ref . [s^l • 

The bottom panel of Fig. [5] shows a situation, where the 
voltage is set to V — 26T. Here, the blocking effect is 
clearly visible since the current is only slightly changing 
in time but has completely different values depending 
on the initial preparation even for times of the order of 
lOr^^. The influence of the phonon shake up process is 
becoming smaller since more phonon modes contribute 
due to the increased voltage window. 
The polaron blocking effect observed in the current 
should also be present in the dot occupancy. In fact, 
within the APTA the two quantities are intimately con- 
nected. In Fig.[n]the dot population for V = 26r is shown 
for diagMC and APTA. For the time accessed here, the 
dot occupation obtained from APTA is only changing 
slightly so that the initial configuration is preserved even 
for times of the order of < > lOF^^. 
The numerical data for the time dependent dot occupa- 
tion of the initially occupied dot are matching with a 
high accuracy the APTA results. On the other hand, 
clear deviations can be seen for the initially empty dot. 
The APTA seems to overestimate in this case the time 



scale for the evolution of the charge. 
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FIG. 6: The dot occupation for the diagMC and APTA are 
compared for V = 26r. While for the dot occupation of 
the initially occupied dot a good agreement is observed, the 
APTA dot occupation of the initially empty dot is overesti- 
mating the blocking effect with respect to the exact diagMC 
results. 

One reason for this deficiency is that PTA underestimates 
the width of the resonances in the spectral density far 
from the Fermi energj*2i. This deficiency is less impor- 
tant for the evaluation of the current as it is determined 
by the resonances within the energy window imposed by 
the electrodes' chemical potentials. An additional source 
for the discrepancy can be related to the finite bandwidth 
which was used for the numerical simulation. In contrast 
to the analytical approach, for a numerical evaluation it 
is necessary to truncate the density of states in the leads 
at some finite energy. Since electron transport far away 
from the Fermi level is important for the time dependent 
dot occupation it can be strongly infiuenced by such a 
finite energy cutoff. This explanation is consistent with 
the findings of Ref. 0] where a strong dependence of the 
time dependent dot occupation on the size of the band- 
width was seen in the Anderson impurity model. 
Further increasing the voltage to = 40r multi-polaron 
processes become more important. This leads to small 
deviation between the APTA results and the diagMC 
data as it can be seen in Fig. [71 Still, the APTA pro- 
vides a qualitative description showing that the transient 
currents neither reach a joint steady state, nor a plateau 
value within the time of several lOF"^. 
Finally a.t V — SOF the deviation between the two cal- 
culations becomes more pronounced. While the diagMC 
results for this bias voltage suggest a convergence towards 
a steady state on shorter time scales, the APTA for the 
initially empty dot exhibits a slower convergence. On the 
other hand, the APTA current for the initially occupied 
dot still reproduces fairly well the numerically exact re- 
sult. 

In any case it should be remarked that the APTA pre- 
dicts in general a slower convergence to the steady state 
than the diagMC results as it can be seen in Figs. |6] and [T] 
This deficiency, which is more pronounced at larger volt- 
ages, can be traced to the already mentioned limitation 
of the PTA spectral density. 
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FIG. 7: The current from diagMC and APTA are compared 
with the same color code and parameters as in Fig.[4]but with 
different voltages, V = 40r (top panel) and V = SOF (bottom 
panel) . 



B. Long time scales 



ferent initial preparations are separated at times much 
larger than several . The polaron blocking effect is 
clearly visible leading to a slowly varying and almost 
constant current for each initial occupation up to times 
t « lOOr^^. For larger times, the polaron blocking is no 
longer the dominant effect, which leads to a charging of 
the dot and the two currents start to converge. Finally 
a joint steady state is reached for much larger times (See 
inset in Fig. [51). 




200 300 



3000 



500 



FIG. 9: Time dependent dot occupation for the same parame- 
ters as in Fig. [S] The charge is blocked for the initially empty 
dot, leading to a small change in the dot occupation up to 
t ~ 100r~^. Inset: The same plot for a larger time scale. 



For the short time scales analyzed in the previous sec- 
tion it was shown that the phonons induce a blocking 
effect which leads to a different time-evolution depend- 
ing on the initial preparation. Further, we have shown 
that for small to intermediate voltages the APTA de- 
scribes correctly the transient behavior of the current for 
the times accessible by the numerically exact diagMC. 
In this subsection we use APTA to analyze the long time 
behavior inaccessible to numerically exact simulations. 




" 

I I I I I I I 

50 100 150 

t [r-1] 

FIG. 8: Time dependent current for long times with the same 
parameters as in Fig. 0] and with V = SF. A plateau value 
of the current for the initially empty dot is observed between 
t ~ 20 — lOOF"'^. Inset: The same plot but for much longer 
times. 

We first focus on the current for the case V = 5T which is 
plotted in Fig. [51 Here, the transient currents from dif- 



The evolution of the dot occupations for the same val- 
ues of the parameters as in Fig. [51 is shown in Fig. [HI 
The charge of the initially occupied dot is close to its 
steady state value so that the transient behavior is not 
pronounced. In contrast, for the initially empty dot the 
occupation is slowly changing until t ~ lOOF"-'^. Then a 
rapid increase is observed and the dot occupation tends 
towards its stationary state value. 

This behavior can be understood from Eqs. (|26l) and (|27p . 
giving the evolution of the charge with time. In fact, since 
F; ^ F, nD(i) can be approximated by 

1 ^ ~ 2 

+ fEf^W(l-^"'''*'0 ' (32) 
/=o 

which explicitly exhibits the fact that only one phonon 
resonance (corresponding to the term in Fj~) lies within 
the voltage window in this small voltage range {V ^ 5F). 
On the other hand, the terms proportional to T^{t) arise 
from the occupied resonances below this window. In the 
case of an initially empty dot, the charge starts to in- 
crease with time dominated by the term in T^it) which 
behaves as ri(l — n-£,{t))^t'^ at short times. The other 
terms in the above equation give contributions propor- 
tional to nD(i)^i^ and are therefore negligible at initial 
times. When time increases an exponential regime is 
reached when the terms in (t) become important. This 
change is rather abrupt and happens at times of the order 
of T'l'{t)t ^ 1, which roughly corresponds to < ~ lOOF""'^ 
for the dominant term I = A. As the voltage increases and 
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more resonances enter in the voltage window, the charge 
of the dot increases more quickly and therefore the transi- 
tion to the exponential behavior occurs at shorter times. 
This leads to a highly non-linear behavior of the charge 
both as a function of time and voltage. 
The currents for V ~ 26r arc showing a similar long 
time behavior (see Fig. fTU]) . After the current oscillations 
from the phonon shake up process die out, the current 
remains almost constant, it only changes slightly in time 
until t « 50r^^ is reached. For larger times a relatively 
rapid charging of the dot causes the currents to finally 
reach their joint steady state in a similar way as in the 
case with V = SF. 




25 50 75 100 

t [r-1] 



FIG. 10: Time dependent current with the same parameters 
as in Fig. |8] but with V — 2GV. The length of the plateau 
is getting smaller, but still it is clearly visible between t ~ 
10 — 50r^^. Inset: The same plot for larger times. 

Increasing the voltage to ^ = 40F the convergence to the 
steady state becomes rather monotonous. Still, the time 
scales involved are much larger than expected for a pure 
electronic system [3] . For even larger voltages this time 
scale is further reduced, however, as commented in the 
previous section, the PTA would not be able to describe 
this large bias regime accurately. 



^2 -V 

^1 ■ — 

r^'^ 

' • ' • ' • ' • ' • ' 

100 200 300 400 500 



FIG. 11: Time dependent current with the same parameters 
as in Fig. [TO] but with V = 40r. The current does not show 
a plateau for small times but the times until a joint steady 
state is reached is still large. 



Finally, we like to make contact between our theoretical 
findings and experiments of single molecular junctions 
by providing a coarse estimate of the set of parameters 
were we expect that such long transients can be found. 
Typical values of F vary between a few fieVs and several 
me Vs. As an example if we set F — meV, the param- 
eters would have the values A = 16 meV, ujq = 8 meV 
and eD = — 10 meV. The applied bias voltages for which 
the long time scales are found would then vary between 
V — 5 meV and V = 40 meV. Accordingly, we would 
obtain transient times of around 50 picoseconds. 



V. CONCLUDING REMARKS 

We have demonstrated the accuracy of the time de- 
pendent APTA method in the strong polaronic regime 
by means of a comparison with the numerically exact di- 
agMC method. A blocking of the current depending on 
the initial occupation for times of the order of several 
F^^ was found in agreement to the results of Ref. [37| . 
Furthermore, a remarkable agreement with these time 
dependent results was found up to moderate voltages. 
We also used this method to explore the long time scales 
which are inaccessible for the exact numerical calcula- 
tions. The polaron blocking effect was shown to be con- 
nected to the narrowing of the side bands in the spectral 
density, determined by F; = Fe~^^ instead of the un- 
renormalized width F. In this way the time scales of 
the system can increase by more than one order of mag- 
nitude in the polaronic regime. Increasing the voltage, 
additional side bands contribute to the electronic trans- 
port which leads to a faster convergence to the steady 
state. 

Finally, some limitations of the method developed in this 
work should be mentioned. Already in the equilibrium 
case the PTA spectral density underestimates the width 
of the side bands far from the Fermi level. In a simi- 
lar way, when a large bias is applied inelastic processes 
not included in the approximation would become impor- 
tant leading to a faster convergence to the steady state. 
Further work along this line would be desirable. 
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